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Abstract 


Artificial luirnorical dissipation is an important issue in large Reynolds number 
computations. In such computations, the artificial dissipation inherent in traditional 
numerical schemes can overwhelm the physical dissipation and yield inaccurate results 
on meshes of practical size. In the present work, the space-time conservation element 
and solution element method is used to construct new and accurate implicit numerical 
schemes such that artificial numerical dissipation will not overwhelm physical dissipa- 
tion. Specifically, these schemes have the property that numerical dissipation vanishes 
when the physical viscosity goes to zero. These new schemes therefore accurately 
model the physical dissipation even when it is extremely small. The new schemes pre- 
sented are two highly accurate implicit solvers for a convection-diffusion equation. The 
two schemes become identical in the pure convection case, and in the pure diffusion 
case. The implicit schemes are applicable over the whole Reynolds number range, firom 
purely diffusive equations to convection-dominated equations with very small viscosity. 

The stability and consistency of the schemes are analysed, and some numerical results 
are presented. It is shown that, in the inviscid case, the new schemes become explicit 
and their amplification factors are identical to those of the Leapfrog scheme. On the 
other hand, in the pure diffusion case, their principal amplification factor becomes the 
amplification factor of the Crank-Nicolson scheme. 

1 Introduction 

The method of space-time conservation element and solution element (the CE/SE method, 
for short) is a new niunericaJ discretization method for solving conservation laws [1-14]. 
It aims to overcome the limitations of the established methods, and is designed to be a 
mathematically simple yet general, accurate and robust method for solving conservation 
laws. The emphasis is on the simulation of the integral forms of the laws, rather than the 
differential forms. This allows for better simulation of regions of rapid change, such as shocks 
and boundEiry layers, where the numerical solution is not smooth. The method has been 
developed on the basis of local and global flux conservation in a space-time domain, in which 
space and time are treated in a unified manner. Derived properties of the conservation laws, 
such as the chaffacteristics and the shock-jump conditions, are not used in the construction 
of the method. For problems in multiple spatial dimensions, the method is genuinely 
multidimensional and is fully compatible with unstructured meshes [2, 4). The CE/SE 
method was first published in this journal in 1995 [Ij. The generalization to multiple spatial 
dimensions appeared in this journal in 1999 [2]. The Introduction section of [1] begins with 
a lengthy discussion that focuses more on physics than on numerics. FVom this discussion, 
readers can understand (i) the considerations that motivate the new method, and (ii) the 
key differences that separate it from the established methods. 

By using a set of design principles that are extracted from the discussion mentioned 
above, several two-level explicit schemes were constructed in [1, 10] to solve (i) the pure 
convection equation 


dufdt -f adujdx = 0 


( 1 . 1 ) 
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and (ii) the convection-diffusion eciuation 

dujdt + a du/dx — fi d\/ dx^ = 0 (1-2) 

where the convection velocity a, and tlie viscosity coefficient fjL (> 0) are constants. These 
schemes were then extended to solve the 1-D time-dependent Euler and Navier-Stokes equa- 
tions of a perfect gas [1, 10]. Moreover, the above 1-D solvers of Eq. (1.1) have been 
generalized to their 2-D counterparts [2, 14, 12]. Because of the inherent simplicity and gen- 
erality of the current method, the above multidimensional generalization is a straightforward 
matter. Also, as a result of the similarity in their designs, each of the above 2-D schemes 
shares with its 1-D version virtually the same fundamental characteristics. 

In this paper we shall describe a simple and innovative approach by which accurate 
implicit time-marching solvers can be constructed using the CE/SE method. A portion of 
this work was presented in [5], and a more general version was recorded in [6]. A striking 
feature of this new treatment is that the modeling of the diffusion-related terms involves 
interpolation between neighboring mesh points while that of the convection-related term 
does not. As a preliminary, first we shall discuss the pros eind cons of explicit and implicit 
schemes. 

For a two-level explicit scheme, the value of a solution at any mesh point has a finite 
domain of dependence at the previous time level. As an example, consider a finite-difference 
solver for Eq. (1.1). Let u", the mesh value of u at any mesh point (j,n) (point P in 
Fig. 1(a)), be determined by Ujl}, and Uj~}. Then the domain of dependence of Uj at 
the (n — l)th time level contains three mesh points. Also one can see that u" is dependent 
only on the initial data given on the line segment AB. 

For an initial- value problem, such as a time-dependent Erder problem or a problem 
involving Eq. (1.1), the solution at any point in space-time also has a finite domain of 
dependence on the initial plane. As a result, explicit schemes could be ideal solvers for such 
a problem if they satisfy the requirement that the physical domain of dependence be a subset 
of the numerical domain of dependence. 

On the other hand, the solution of an initial-value/boundary-value problem at any point 
in space-time is dependent on the initial data and the boundary data up to the time of the 
point under consideration. As an example, consider a problem involving Eq. (1-2). As shown 
in Fig. 1(b), the solution at point P is dependent on the initial/bormdary data given on EC, 
CD, and DF where E, P, and F are at the same time level. Let this problem be solved 
using the explicit scheme that was explained using Fig. 1(a). Let P also be a mesh point 
(j,n). Then uJ is dependent only on the initial/boundary data given on AC, CD and DB. 
It is completely independent of those data, given on AE and BF. Contrarily, if the same 
problem is solved using an implicit scheme, then u” is dependent on the initial/boundary 
data given on EC, CD, and DF. In other words, the numerical domain of dependence of 
the implicit scheme is consistent with the physical domain of dependence of the problem 
under consideration. 

Two observations can be made as a result of the above discussions. 

(a) Generally an explicit scheme is not an ideal solver for an initial-value/boundary-value 
problem. Because a time-dependent Navier-Stokes problem is such a problem, the 
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iibove arguinoat. implies that an explicit scheme cannot be used to solve? a time- 
dependent Navi(?r-Stokes problem except for the special circumstaiure in which errors 
caiLsed by neglecting certain initial/boundary data (such tis those given on AE and 
BF \n Fig. 1(b)) are relatively small. The factors that help achieve the above special 
circumstance include: (i) a small time-step size to spatial-mesh interval ratio, (ii) a 
small time rate of change of boundary data, and (iii) a small contribution of the viscous 
terms in the Navier-Stokes equations relative to that of the inertial terms. Note that 
condition (iii) may be met by a high-Reynolds-number flow. 

(b) Generally an implicit scheme is not an ideal solver for an initial- value problem. This is 
because the domain of dependence of the former is greater than that of the latter and, 
as a result, an implicit solution tends to be contaminated by extraneous information. 


2 Preliminaries 

Let Eq. (1.1) be in a dimensionless form. Let xi = x and X 2 = t\)e considered as the 
coordinates of a two-dimensional Euclidean space By using Gauss’ divergence theorem 
in the space-time E^, it can be shown that Eq. (1.1) is the differential form of the integral 
conservation law 

<f h-ds = 0 (2.1) 

Js(V) 

Here (i) ^(V) is the boundary of an arbitrary space-time region V in E 2 y (ii) h — (ou, u) 
is a current density vector in smd (iii) ds = daH with da and n, respectively, being 
the area and the outward unit normal of a surface element on S{V). Note that (i) h • da 
is the space-time flux of h leaving the region V through the surface element ds, and (ii) all 
mathematical operations can be carried out as though were an ordinary two-dimensional 
Euclidean space. 

In the following, we shall briefly review the inviscid version of the explicit a-p, scheme 
[1, 10]. Let Hi denote the set of mesh points (j,n) in (dots in Fig. 2) where n = 
0, ±1, ±2, ±3, . . . , and, for each n, j = n, n ± 2, n ± 4, . . . . There is a solution element (SE) 
associated with each (j,n) € fii- I^t the solution element SE(J,n) be the space-time region 
boimded by the dashed curve depicted in Fig. 3. It includes a horizontal line segment, a 
vertical line segment, amd their immediate neighborhood. 

For any (x, t) G SE(j, n), u{x, t), and h(x, t), respectively, are approximated by iz*(x, t ; j, n) 
and h*{x,t;j,n), which we shall define shortly. Let 

ti*(x, t]j, n) = + (u,)?(x - Xj) -b {ut)j{t - r) (2.2) 

where (i) (u,)^, and (ut)y are constants in SE(j,n), and (ii) (xj,t") are the coordinates 

of the mesh point {j,n). Note that u", (u^)", and (ut)" can be interpreted as the numerical 
analogues of the vaJues of u, dufdx, and du/dt at (xj,t'‘), respectively. 

We shall require that u = u*{x,t]j,n) satisfy Eq. (1.1) within SE{j,n), i.e., 

(nO” = -a(ux)" (2.3) 
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Coinhining Eqs. (2.2) aiid (2.3), ono has 

u*{x, t -J, n) = Uj + {u^yj [(x - xj) -a{t- f")] (x, t) € SE(j, n) (2.4) 

Because h = {au, u), we define 

h*{x,t-,j,n) = {au*{x,t]j,n), u{x,t]j,n)) (2.5) 

Let E 2 be divided into nonoverlapping rectangular regions (see Fig. 2) referred to as 
conservation elements (CEs). As depicted in Fig. 4, the CE with its top- right (top-left) 
vertex being the mesh point (j,n) 6 is denoted by CE_(j,n) (CE+(j, n)). Obviously the 
boundary of CE_(j,n) (CE+(j', n)) is formed by subsets of SE{j,n) and SE(_; — l,n — 1) 
(SE(j 4- l,n — 1)). The current approximation of Eq. (2.1) is 

F±{j,n)= (f h*-ds = 0 (2.6) 

for all (j,n) E Qj. In other words, the total flux leaving the boundary of any CE is zero. 

Because (i) The CEs associated with Oi can fill any space-time region, and (ii) the 
surface integration across any interface separating two neighboring CEs is evaluated using 
the information from a single SE, the local conservation condition Eq. (2.6) leads to a global 
conservation relation, i.e., the total flux leaving the boundary of any space-time region that 
is the union of any combination of CEs will also vanish. 

With the aid of Eqs. (2.4)-(2.6), it can be shown that 

F^U, n)/Ax = ±(1 - i-’) [«)” + («:)-■] + (1 TO K - «Kl‘) (2-7) 

where u = aAt/Ax is the Courant number, and (itj)" = {Ax/2){ux)j. Note that here Ax 
and At, respectively, represent the same mesh interval and time-step size which were denoted 
by Ax/2 and At/2 in [1, 10]. Using Eqs. (2.6) and (2.7), u" and (uj)”, which are considered 
as independent unknowns at the mesh point {j,n), can be solved for in terms of u^f.\ and 
if 1 — 0. It can be shown that, for all (j,n) € fli, 

q{j,n) = Q+q{j -l,n-l)-^Q.q{j-\-l,n-l) (1 - ^ 0) (2.8) 

Here (i) q(j,n) is the column matrix formed by and («J)", and (ii) 

Q± = (1/2) (^ ) (2.9) 

Eq. (2.8) defines a marching scheme. Because this scheme is the special case of the a-p 
scheme when p = 0, hereafter it will be referred to as the a scheme. It is the only two- 
level explicit solver of Eq. (1.1) known to the authors to be neutrally stable, i.e., free from 
numerical dissipation. 

In the above construction of the a scheme, we use the SEs and CEs of the mesh points 
marked by dots in Fig. 2. A similar construction can be performed by using the mesh points 
marked by triangles in Fig. 2. Let ^2 denote the set of mesh points {j, n) in E^ (triangles 


4 




in Fig. 2) whoro n = 0, ±1, ±2, ±3, . . . , and, for each n, j = n. ± 1, /?. ± 3, n ±5, — Let 
the; SEs and CE.s c^f Q.) ix; defiiK'd by using Figs 3 and 4 with dots replaced by triangles. 
Obviously (i) the CE.s of Q -2 also fill any space-time region, and (ii) the a sc:hernc can also 
be constructed using the SEs and CEs of ^ 2 . This new scheme is defined by E((. (2.8) with 

(j i ^ ) ^ ^2- 

Before we proceed further, some intricate points related to the above constructions will 
be clarified with the following remarks: 

(a) SE(j, n) may intersect SE(j',n') if {j,n) € and 6 ^ 2 - 

(b) Let {j,n) € Qi. Then (i) (j -j- l,n) € Ct 2 , smd (ii) CE+(j,n) and CE_(j -1- l,n) represent 

the same rectangle in E 2 . However, because the function h* used in the evaluation of 
F+{j,n) is tied to a pair of SEs associated with fii, while that used in the evaluation 
of F_(j -1- l,n) is tied to another pair of SEs associated with fl 2 ) = 0 and 

F_(j d- l,n) = 0 represent two completely independent flux conservation conditions. 

To prepare for the development of the implicit solver to be described in the next section, 
we shall combine the above two independent schemes into a single scheme. The new scheme, 
referred to as the dual a scheme, is defined by Eq. (2.8) with {j, n) eQ where = f2i U f22- 
Obviously, a solution of the dual a scheme is formed by two decoupled solutions with each 
being associated with a mesh that is staggered in time. Several classical schemes also have 
this property. Among them are the Leapfrog, the DuFort-fVankel, and the Lax schemes [15]. 

3 The implicit schemes 

An implicit solver for Eq. (1.2), referred to as the a-p,{Il) scheme, will be discussed first 
in this section. Here stands for “implicit”, and “1” is the identification number. This 
discussion will be followed by a brief discussion of another implicit scheme, termed the a- 
scheme. The implicit schemes are constructed to meet two requirements given in the 
following discussion: 

(a) With a few exceptions, numerical dissipation generally appears in a numerical solution 

of a time-marching problem. In other words, the numerical solution dissipates faster 
than the corresponding physical solution. For a nearly inviscid problem, e.g., flow at a 
large Reynolds number, this could be a serious difficulty because numerical dissipation 
may overwhelm physical dissipation and cause a complete distortion of solutions. To 
avoid such a difficulty, the model solver is required to have the property that the 
numerical dissipation will approach zero as the physical dissipation approaches zero. 

(b) The convection term and the diffusion term in Eq. (1.2) involve the spatial derivatives of 

first order and second order, respectively. Thus, in a spatial region where a solution is 
very smooth, the diffusion term is negligible compared with the convection term. As a 
result, the effective physical domain of dependence is more or less dictated by Eq. (1.1). 
To prevent excessive contamination of the solution by extraneous information, the 
implicit solver mil be required to become an explicit solver in the limiting case in which 
the diffusion term vanishes. 
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Because of the rcquireiiKuits set forth in (a) and (i)), the implicit solver will be constructcxl 
such that it reduccss to th(; dual n schenu! if ji = 0. Tlui foriiK^r diffiTS from the latter only in 
th<! extra modeling involving the diffusion-related terms. Note that the presence of viscosity 
is felt through (i) the diffusion term -fxd'^ufdx^ in Eq. (1.2), and (ii) the spatial diffusion 
flux component —fxdu/dx in the flux vector 

h = {au — ndu/dx, u) (3-1) 


Note that Eq. (1.2) is the differential form of Eq. (2.1) if h is defined according to Eq. (3.1). 
Also, in this paper, any term in Eq. (1.2) or on the right side of Eq. (3.1) is considered to 
be convection-related if it is not diffusion-related. 

To construct the a-/z(/l) scheme, consider a finite portion of the mesh depicted in Fig. 2 
(J > 4). We assume that (i) u = uj(x) at f = 0, (ii) u = ui{t) at a: = 0, and (iii) u = Ufi{t) 

at X = JAx, where u/(x), U£,(t), and are the given initial data, left-boundary data 

and right-boundary data, respectively. Moreover, for the current case, (i) and 0.2 are 
restricted by the conditions n > 0 and J > j > 0, (ii) CE±(j, n) are not defined if n = 0, 

(iii) CE_(j,n) is not defined if j = 0, and (iv) CE+(j,n) is not defined if j = J. Items 

(iii) and (iv) imply that only one conservation condition is associated with a boundary mesh 
point. Obviously, the definition of SE(j, n) also needs to be appropriately modified if j = 0, 
or j = J, or n = 0. 

Eq. (2.2) will still be assumed. We also assume that, for n = 0, 1, 2, ... , 


y^=ULin yri=nRin («o}=«n(n (3.2) 

where Ui(t) dui{t)/dt and UR{t) = duR(t)/dt. Thus, only one unknown, i.e., (ui)”, and 
one conservation condition are associated with a boundary mesh point {j, n). 

Furthermore, for an interior mesh point (j,n), we replace Eq. (2.3) with 


(f,)" = -a(«.)J + ^ [(u.)",., - («.);_,] J>i>0 (3.3) 

Eq. (3.3) is the numerical analogue of the differential condition Eq. (1.2). Eqs. (2.2) and 
(3.3) imply that, for J > j > 0 and (x, t) 6 SE(j,n), 

u*{x, t;j, n) = u] + (Ux)j (x ~ Xj) + [M]+i ~ («x)”_i] - a(«x)7 } (^ “ (3-4) 

Next, as a result of Eq. (3.1), Eq. (2.5) is replaced by 

h*{x,t;j,n) = {au*{x,t\j,n) - iml{x,t\j,n), u*(x,t\j,n)) (3.5) 

Here, for any (x,t) €SE(j,n), 







/2, ift>t”; 
/2, ift<t". 


(3.6) 


Consider any point (x, t) on the line segment joining the mesh points (j, n) and (j, n - 1). 
Then {x,i) belongs to both SE{j,n) and SE(j, n — 1). According to Eq. (3.6), for this point 
(x, t), 


ul{x,t]j,n) = ul(x,t;j,n - 1) = [M] + M] /2 


(3.7) 
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Thus, the same numerical diffusion flux component is assigned to the point {x,t) regardless 
of whether it is considered ns a point in SE{jca) or a point in SE{j,n — 1), in tlieir r(?gion 
of overlap. 

With the above modifications, the a-p[1 1) scheme is defined by assuming Eq. (2.6). Note 
that: 


(a) At a mesh point e the diffusion-related terms in Eqs. (1.2) and (3.1) are modeled 

using inteTT)olations that may involve the numerical values of the mesh points G fl 2 (^i)- 
This contrasts sharply with the modeling of the convection-related terms which uses no 
interpolation. 

(b) In the dual a scheme, the two sets of numerical variables associated with fli and ^2 

are completely decoupled from each other. Contraxily, they are "glued” together in the 
a-^(Il) scheme through the interpolations referred to in (a). 


To proceed, let a = fj,At/{Axff, and {uf)f (At/2)(ut)”. Also let 


III 

(1 - - (1 - - a)(Wx)j+i 



jj/y 

- (n = 0, 1,2,. . . = 0, 1,2, . . . , J - 2) 

(3.8) 

1 

III 

i - a(“?)3-i - (1 - - "W); (n = 0, 1,2, . . .) 

(3.9) 

{S-)f = 

(1 + + (| + l)a(0; + (1 - 1/^ - a)K);_, 



i/rv 

(n = 0, 1, 2, . , . ; j = 2, 3, 4, . . . , J) 

(3.10) 

= (1 + t'K 

+ a(uj )f + (1 - q){uJ)S + i/(«+)S (n = 0. 1,2, . . . ) 

(3.11) 


By using Eqs. (2.2), (3.2), (3.4)-(3.6), and (3.8)-(3.11), Eq. (2.6) implies that 

(1 T l^)u^ ± (1 - + a)«)" + (j T l)a«)"±i - 

= (S±)r' = = 


(3.12) 


(1 + a)(uj); - = (5Ur‘ - (1 - "K - ■'W)S (r. = l,2,...) (3.13) 

«KU"-,-(i+“){“J)3 = (s-)r‘-(i+>')«;+KO; (n=i,2,...) (3.u) 

Eqs. (3.12)-(3.14) define the a-p.{I\) scheme. Note that Eq. (3.12) represents a pair of 
equations for each (j, n). 
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This section is roncliidcd with a discussion of another implicit schorno, referred to t.\s the 
a-n{12) schimie. It difhns from the a-/i(/l) scheme only in one n'spect, i.e., Eq. (3.3) is 
replaced by 

= + J>J>0 (3.15) 

As a result, Eq. ( 3 . 4 ) must be modified accordingly. Note that (i) because the a-/i(/l) and 
a-/x(/ 2 ) schemes differ only in the representation of the viscous flux in Eqs. (3.3) and (3.15), 
the two schemes are identical if /i = 0 , and (ii) because the differing expressions for (ut)" in 
Eqs. (3.3) and (3.15) are applied only in the evaluation of the flux of au* in the conservation 
equations, the two schemes are identical if a = 0. Let 

(fl+)J = (1 - >- + - “(“i")" - (1 - - a)(“J)”+i 

-^(“"+< 2 ) (n = 0,l,2....;j=0,l,2,...,J-2) (3.16) 


(1 - ^)«5 _ a(uj)5_i - (1 - «)«)!} - (n = 0, 1,2, . . . ) (3.17) 


(i?_)” = {l + u- ua)u]_, + ««)” + (1 - i/" - a)«)”_i 

+l^{u] + u]_ 2 ) (n = 0,l,2,... =2,3,4,... ,J) (3.18) 

(i?_)? '= (1 + + a{utr, + (1 - a)«)J + u{ut)^ (n = 0, 1, 2, . . . ) (3.19) 

Then in the a-fi{12) scheme, Eqs. (3.12)-(3.14) are replaced by 

(1 ip uct)vr^ + a)(n+)” T a(«^)"±i ± ^ (u^+i + u"_i) 

= (i?±),"-' (n=l,2,3,...;j = l,2,...,J-l) (3.20) 

(1 + a)«)o - a(Or = {R+)V - (1 - (n = 1, 2, . . . ) (3.21) 

««)”-! - (1 + a){ut)j = {R-)T^ - (1 + i'Kj + HO" (n = 1,2 , . . . ) (3.22) 

Eqs. (3.20)-(3.22) define the a-^(/2) scheme. Note that Eq. (3.20) represents a pair of 
equations for each (j, n). 


4 Solution Procedure 

We first discuss the solution procedure of the a-/x(/l) scheme. Let 1 — 7 ^ 0. Then 

Eq. (3.12) is equivalent to the pair of equations 

+ 2(1 c )«)7 - «(<),"« = ( 1 + i-)(sur‘ - (> - •')(s-)r' (-‘.i) 
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and 


“ 5 {(s+)r‘ + (s-)r' + “ [«)"+! - 

where j = 1,2,3,.. . , 7 - 1. In the following discussion, Eqs. (3.12)-(3.14) will be replaced 
by Eqs. (3.13), (3.14), (4.1) and (4.2). 

Let the marching variables at the (n — l)th time level be given. With the aid of Eq. (3.2), 
the expressions on the right sides of Eqs. (3.13), (3.14) and (4.1) can be considered as given 
source terms. Thus these equations form a tridiagonal system of J + 1 equations for the 
J + 1 unknowns (u+)^, ; = 0, 1, 2, . . . , J. It can be shown [6] that the coefficient matrix 
associated with the tridiagonal system is strictly diagonally dominant in rows and columns, 
i.e., stability of the Thomas algorithm ([15], p.99) for solving the system is assured, if 

i/^ < 1 and a > 0 (4.3) 

Upon obtaining j = 0, 1,2, . . . , J by solving the tridiagonal system, the other un- 

knowns j = 1, 2, . . . , J - 1, can be obtained using Eq. (4.2). The time-marching stability 
of the a-fjL(Il) scheme, as well as that of the o-/i(/2) scheme, is discussed in the next section. 

Finally, we briefly discuss the solution procedure of the a-fi{12) scheme. If the marching 
variables at the (n - l)th time level and the boundary conditions Eq. (3.2) are given, the 
expressions on the right sides of Eqs. (3.20)-(3.22) can be considered as given source terms. 
Thus these equations form a system of equations for the u" and the {u^)j at the nth time 
level. Unlike the case with the a-/x(/l) scheme, prior elimination of the u" at each mesh 
point to obtain a system for only the (u^)" is not possible, as the conservation equations 
(3.20) are implicit in both and (u+)”. Hence, the resulting system for the and the 
(u+)" must be solved by a block version of the Thomas algorithm. This makes the a-^(/l) 
scheme computationally more efficient than the a-n{12) scheme. 


5 Stability Analysis 



\ r 

n) = , 

L Wx)j J 

(5.1) 

qL‘| 

/ 0 —va/2 \ 

0 (i//2-bl)a ) 

(5.2) 

,(1) ( l-v l-u^ + a \ 

° — (1 — i/^ + a)y 

(5.3) 

of •n 

(0 (i//2-l)a'\ 
1 0 —i/aj2 J 

(5.4) 
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^ t 


't/ 


0 0 \ 
0 —ua/2 J 


(5.5) 


q(2) */ 


0 0 
1 + 1 / 1 — i/^ — a 


(5.6) 


_(2) de/ / 0 (t^/2 1) a \ f5 7l 

^0 “ V 0 (j//2 + 1)q y ^ 

-( 17 “-“)) (5.8) 


gr = 


Then Eq. (3.12) may be expressed as 


(2) d^f ( 0 -uotl2 \ 

Vo 0 y 


Q?~^{j + hn)=Y “ 1)- 


(5.9) 


(5.10) 


We examine the evolution of a single Fourier mode of the error by making the substitution 
n) = q*{n, 6) e*-'^ (i = -tt <0 <ir) (5.11) 


in Eq. (5.10), where (i) q*{n,$) is a 2 x 1 column matrix, and (ii) 9, —tt < 6 < tt, is the 
phase angle variation per unit Ax of a single Fourier mode. Performing this substitution, 
we obtain 


where 


and 


Q^'\v,a,e) q’(n,e) =Q‘-'‘\v,a,9) q'{n- l,$), 




1 - 1 / 

1 + 1 / 


i=-l 

1 - 1/2 + a - if e-‘® + Q (| - 1) e*'' 

- (1 - 1/2 + a) - f + a (f + 1) e"*® 


(5.12) 


(5.13) 


/=-2 


(1 — z/)e'® 

(1 + u)e~'^ 


- (1 - 1/2 - o) e'® - if e2‘<» + a(^-l) ' 

(1 - i/2 - a) c-^ - if e-2*i' + a (^ + 1) ’ 


(5.14) 
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'= det. 

= -2 [(1 - +a(l -COS0)] , (5.15) 

it is evident that 7 ^ 0 for — tt < 0 < tt if the conditions Eq. (4.3) are satisfied. Assuming 
Eq. (4.3), then exists and we can multiply Eq. (5.12) from the left by 

As a result, the amplification matrix is 

G(u,a,e) = [Q^^\i^,a,0)]~'Q^'^\u,a,e). (5.16) 

The amplification factors of the 1) scheme are the eigenvalues of G. It is straightforward 

to show that any eigenvalue X{u,a,6) of G{u,a,9) satisfies the condition 


o 

II 

1 

(5.17) 

Using the definitions Eqs. (5.13) and (5.14), Elq. (5.17) is equivalent to 


AX^ + BX + C = 0, 

(5.18) 

where 


A = (1 - -1- O' (1 - COS0) 

(5.19) 

5 = 2a (l - cos 0 - r/^ sin^ 9) + 2iv (l - i/^) sin9 

(5.20) 

and 


C = - (1 - u^)+a{l-cos9). 

(5.21) 

Thus, the amplification factors A± are given by 

, -B ± - 4AC 

2A • 

(5.22) 


Note that there are two amplification factors rather than one, because there are two un- 
knowns at each mesh point. Also note that (i) A± — » ±1 and 0 — > 0 if 1 — 7 ^ 0, and 

(ii) the analytical amplification factor of a plane- wave solution to Eq. (1.2) approaches 1 
as 0 ► 0. Therefore, A+ and A_ are referred to as the principal and spurious amplification 

factors, respectively. Numerical evaluations of Aj. have shown that the scheme is 

stable provided that the conditions Eq. (4.3) are satisfied. 

Note that many other implicit solvers are unconditionally stable. However, the price paid 
for this “desirable” property usually is excessive numerical dissipation. Moreover, the use 
of a time-step size that is greater than that allowed by Eq. (4.3) generally results in a less 
accurate time-dependent solution. Thus we do not consider the more restrictive stability 
condition Eq. (4.3) to be a disadvantage of the a-/i(/l) scheme. 
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When ^ = 0, the a-ii{Il) scheme reduces by design to tlu; dual a scheme. Also, Eq. 
(5.22) implies that 


A± = -ii/sin9 ± V 1 - if /i = 0 


(5.23) 


The amplification factors given in Eq. (5.23) are exactly identical to those of the classical 
Leapfrog scheme ([15], p.lOO). On the other hand, Eq. (5.22) implies that 


1 - a (1 — COS0) 

1 + a (1 - cos 6) 


and A_ = — 1 


if a = 0 


(5.24) 


It is found that in this case A+ is identical to the amplification factor of the Crank-Nicolson 
scheme ([15], p. 112). 

The fact that the amplification factors of the a-fi{Il) scheme are related to those of 
two celebrated classical schemes is only one among a string of similar coincidences. Other 
coincidences are summarized in the following remarks [1, 6, 10]: 


(a) By using exactly the same procedure by which the dual a scheme is formed from two 

decoupled a schemes, the dual a-p scheme and the dual a-e scheme can be constructed 
using the a-p scheme [1] and the a-e scheme [1], respectively. It can be shown that the 
amplification factors of the dual a-p scheme reduce to those of the Leapfrog scheme 
when p = 0 and to those of the DuFort-Prankel scheme ([15], p.ll4) when a = 0. 
Also, it can be shown that the amplification factors of the dual a-e scheme (i) reduce 
to those of the Leapfrog scheme when e = 0, and (ii) become the same function of z/ 
and 9 when e = 1 and this function is identical to the amplification factor of the Lax 
scheme. 

(b) It is explained in [1] that each of the Leapfrog, DuFort-Prankel and Lax schemes is 

formed by two decoupled schemes (each of the decoupled schemes is defined on a 
staggered space-time mesh). It is also shown in [1] that the amplification factors of 
the decoupled Leapfrog, DuFort-Prankel and Lax schemes are related to those of the 
a, a-p and a-e schemes in exactly the same ways that the amplification factors of the 
Leapfirog, DuFort-FVankel and Lax schemes are related to those of the dual a, a-p and 
a-e schemes 


In [6], the stability of the a-p{12) scheme is analysed, with the same procedure as was 
\jsed for the a-p{I\') scheme. It is shown there that the a-pi^lTj scheme is stable if the 
conditions Eq. (4.3) are satisfied. As pointed out in Section 3, when a = 0, the a-p{Il) 
scheme and the a-p(I2) scheme become identical. Hence, the a-p{12) scheme shares with 
the a-p{I\) scheme the property that when o = 0, its principal amplification factor becomes 
identical to the amplification factor of the Crank-Nicolson scheme. Since the a-p{12) scheme 
also reduces to the dual a scheme when = 0, its amplification factors also reduce to those 
of the Leapfrog scheme in this case. 
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6 Consistency and Truncation Error 

In this section, w(! discuss the consistency and truncation error of the scheme. For 

a similar analysis of the a-/i(/2) scheme, the reader is referred to [6]. As a preliminary, this 
section will begin with a discussion of a critical concept. 

First, note that in a typical numerical scheme, a physical variable is associated with 
a single numerical variable. Thus, a system of two coupled physical equations involving 
two independent physical solution variables, at each mesh point generally is modeled by 
a system of two coupled discrete equations involving two independent numerical variables. 
Also, one would expect that the two coupled discrete equations are consistent with the two 
coupled physical equations. Thus, in general, one would not expect that two coupled discrete 
equations be consistent with only a single PDE. 

The scheme is nontraditional in one key respect. Even though it is introduced 

to model a single PDE (i.e., Eq. (1.2)) with a single dependent variable u, at each interior 
mesh point it is formed by two coupled discrete equations (see Eq. (3.12)) involving two 
independent numerical variables u” and («i)” (note: (uj )” = Ax/2 (ux)"). 

The numerical variables u" and (ui)" could be “interpreted” as the numerical analogues 
of u and du/dx, respectively. However, it should be understood that this interpretation is 
not exact in nature and that it certainly does not invalidate the fact that uf and (ui)" are 
independent numerical variables. As a result, one would expect that the two equations given 
in Eq. (3.12) be consistent with a system of two PDEs with one of them being Eq. (1.2). 

The interior equations of the a-p{Il) scheme, Eq. (3.12), represent the results of the 
evaluation of Eq. (2.6) with the use of Eq. (3.5). Consider the equations obtained by 
evaluating 


F+{3,n) =0 


( 6 . 1 ) 


and 


F_(i + l,n) = 0 (6.2) 

These are two distinct conservation equations, although the conservation elements involved 
occupy the same physical region. The a-p{I\) scheme can obviously alternatively be de- 
scribed by specifying these two equations. Linear combinations of these equations axe next 
examined, in order to investigate the consistency of the scheme. Adding the equations 
symbolized by Eq. (6.1) and Eq. (6.2), and dividing the result by 2 Ax At, results in 

[FD£l(u,u,)|";|=0 (6.3) 

with FDEl as defined below. Also, subtracting Eq. (6.2) from Eq. (6.1), and dividing the 
result by 2(1 — iP) Ax^, results in 


[FDE2{u,Ux)r.;l = 0 


(6.4) 
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Hero, FDEl and FDE2 are defined as 
|FDEl(u, .))]";! 'tf 


^ Ki - - >*r'i + ^ K. - - “r‘l 


2A^ 


K+i - + <+■.* - ^r'l 


2 Ax 
Ax 

TKi 

VfJ, 


va 

T 


[''i+2 ~ - ^j+2 + 


8 Ax 

iEDE2(t,,ei;;| =' 


(6.5) 


[-JV2 - + »;vi - - >'r'] 


l/Q 


4(1 — 1^2) 
1 


2 (1-1/2) 
V 


b^l+i '“i+i ] 

[u"+i - u]+l + u] - 


( 6 . 6 ) 

2 (1-1/2) Ax . 

Note that, to emphasize the fact that and (ux)" are two independent marching vari- 
ables, a new symbol v is introduced in Eqs. (6.5) and (6.6). 

Let w(x, t) and i;(x, t) be smooth functions. Furthermore, let (i) 


, , def du du d^u 

PEEl(u) 

PDE2{u, v) v-^ 

^ def - 


(6.7) 

( 6 . 8 ) 


and (ii) ii" u (j Ax, nAt) and v v (j Ax, nAt). Let [FDEl(u, v)]”+| and [F DE2(u, u)] 

be considered as the discrete approximations of [PDEl(u)]" ^ and [PDE2(u, v)] i, respec 




j+i 


tively. Then the errors ERl and ER2 in these approximations may be defined by 

[Fi?i];;| [FDFi(u,{5)];;| - [fdfi(u)];;| 

n-i def 


[FF2];;| =' [FFF2(u,5)];;| - lPDE2(u,v)]^_^l 
With the aid of Taylor expansions, it may be shown that 


(6.9) 

(6.10) 


= -4(”- 


d^u Ax2 
+u-;r^“r': !■ ^ 


du\ d^u Ax2 d^uAt^ 
dx'^dt s' " dF 24 
AF d^v Ax2 


+':r 


dx^ 24 
1 d^v 


dx dt'^ 8 


- /2 


dx^ 24 
d^v At^ 


- P 


d^v At- 
dxdi"^ 8 


8 5x dt 


[a^Ai= - Ax^] - + O (A») (6.11) 
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It. may similarly lx; sIkwii tluit 


[ER2] 


j + -2 
u At 


(1 — i/'^) Ax 
d^v Ax^ 
^ax2 8 
ua Ax 
'^12(1 - 1/2) 
1 

~24(1 - //2) 
At 


du du 
dt ^ ^ dx 
d'^v AE 


fJ- 


d-u 

3x2 


+ 


2uaAx 0 
1 — i/2 dx 


( . dn 

V 


3t2 8 

'^V . 2 o 

*^^3x3t2 


Ar 


S^’-sS."’ 


24 (1 - 1/2) Ax 

. fA^l O (A^) O (A^) 

^(1-1/2)^^^ (l-i/2) ^ ^ (l-i/2)Ax ^ ^ 


+ O (A^) 
u At 


( 6 . 12 ) 


In Eqs. (6.11) and (6.12), all derivatives are evaluated at {{j + 1/2) Ax, (n — 1/2) At). Each 
symbol O (A^) represents an infinite sum of terms, in which derivatives of u or v, and the 
quantities a, n and Ax'At'” occur only as factors in the nmnerator of each term, with 
l,m>0 and / + m > 3. 

Let u and ii be a solution of the system of PDEs PDEl{u) = 0 and P DE2{u, v) = 0. 
Then [FDEl{u,v)f^~^ and [FDE2(ti,{i)]”~| are by definition the truncation errors of the 

discrete equations (6.3) and (6.4). Because [PDEl{u)]^^l = 0 and \PDE2{u,v)]^^^ = 0, 

Eqs. (6.9) and (6.10) imply that [FDEl{u,v)\^^^ = [ER\\^^l and \FDE2{u,v)].^\ = 

[Ei?2]""|: Also, because [FI?El(u)]"'| = 0 and [PDE2(ii, u)]”"| = 0, the first term 

on the right hand side of Eq. (6.11) and the first two terms on the right hand side of 
(6.12) vanish. Assume ^ 1. Let the rule of mesh refinement be such that ^ remains 
bounded as Ax — ♦ 0 and At — » 0. This implies that v and a Ax also remain bounded. 
Examination of ER\ and ER2 then shows that (i) the discrete equations (6.3) and (6.4) of 
the a-/i(/l) scheme are consistent with the adv^tipn-diffusion equation PDE\{u) = 0 and 
with PDE2{u,v) = 0, and (ii) the scheme is second-order accurate in space and time. 

In [6], a similar consistency analysis of the a-fi{12) scheme is performed. The conclusion 
drawn there is that if the rule of mesh refinement is restricted as described in the previ- 
ous paragraph, the a-/z(/2) scheme is also consistent with the advection-diffusion equation 
PDE\{u) = 0 and with PDE2{u,v) — 0, and the scheme is also second-order accurate in 
space and time. 


7 Numerical Results 

Three test problems will be used to evaluate the accuracy of the a-/i (71) scheme. It is shown 
in [6] that the a-/z(/2) scheme yields similar results. 


15 


7.1 Decaying Traveling Sine Wave 

In tlio first probl(!in, wo consider a special case of Eq. (1.2) (a = 1 and p = 0.01) with 
0 < a; < 1 and t > 0. The initial and boundary condition functions u/(x), ui,{t) and un{t) 
are defined such that they are consistent with a special solution to Eq. (1.2), i.e., 


u = Ue{x,t) '= exp(-47T^/it) sin [27r(z - at)\ 


(7.1) 


Let Uex(a:, t) '= due{x, t)jdx. At any time t = t”, let 


j-\ 


L,{U) '^ 1 :' 


( J — 1) exp(— 47 tV^") ^ 




(7.2) 


1-1 (Wx) 


1 

(J + l)27rexp(— 47r2/it”) 


j 


j=0 


(7.3) 


Li{u) and Li{u^) are two error norms (per mesh point) which are normalized by the decay 
factors of Ue{x,t) and u^^ix, t), respectively. 

Let J = 80 (i.e., Ax = 1/80) and u = 0.8 (i.e., At = 0.01). Then a numerical computa- 
tion yields Li{u) = 0.7961 x 10“^ at t = 4 (i.e., n = 400). Through numerical experiments, 
it has been shown that both Li{u) and Li(ui) at a given time t are reduced by a factor of 4 
if both Ax and At axe reduced by half, confirming the second-order accurate nature of the 
scheme. 

Fig. 5 shows a comparison of the computed solution at t = 4 with the exact solution. It 
also shows the error Ue{xj,t^) — u" scaled by the peak magnitude of the exact solution at 
that time level. It is seen that the maximum error is less than 0.3% of the peak magnitude. 
The peak magnitude is seen to be just over 0.2. 

Fig. 6 shows a comparison of the errors in the solutions obtained with the a-/z(/l) scheme 
and the implicit MacCormack scheme ([15]). Both schemes were applied with the same 
parametere and mesh as described above. It is seen that the current scheme is considerably 
more accurate than the implicit MacCormack scheme. Refinement of the grid keeping u 
constant would actually make the comparison even more favorable. This is because the 
implicit MacCormack scheme is second-order accurate in space and time only if a is held 
constant and u —* 0 when refining the grid. 


7.2 Pure Diffusion 

We consider a special case of the convection-diffusion equation with a = 0 and = 1, in 
the domain 0 < x < 1 and t>0. The initial/boundary conditions completing the problem 
specification axe (i) u{0,t) = u(l,t) = 0 for f > 0, (ii) u(x,0) = 2x for 0 < x < 0.5, and 
(iii) u(x,0) = 2(1 — x) for 0.5 < x < 1. The solution u{x,t) exhibits the diffusive decay 
of the initial sawtooth shape. An exact series solution is available, see for e.g. p.l5 of [16]. 
For the CE/SE computation, uniform mesh intervals Ax = 0.02 and At = 0.005 axe used. 
Fig. 7 shows the time-slice at t = 0.05, comparing numerical and exact solutions, and also 
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showing the error sealed with the peak exact value at that time level. The maxirnum error 
iiiaguitude is sefai to be about 0.5% of the peak solution value. At /. = 1 (not shown), when 
the peak solution value has dwindled to about 4 x 10“'^, the maximum error magnitude is 
about 0.15% of the peak solution value. 

The same problem on the same mesh was solved with the Crank-Nicolson scheme, which is 
probably the best traditional scheme for the parabolic pure diffusion equation. Fig. 8 shows 
a comparison between the current scheme and the Crank-Nicolson scheme, at t = 0.05. 
The schemes are seen to be of similar accuracy, as is to be expected from the fact that, 
when a = 0, the principal amplification factor of the a-/i(71) scheme is identical to the 
amplification factor of the Crank-Nicolson scheme. 

7.3 Steady State Boundciry Layer 

We next consider the problem defined for the convection-diflfusion equation in the domain 
0 < X < 1 and f > 0 by the conditions (i) u(0,t) = 0 for t >0, (ii) u{l,t) = 1 for t > 0, and 
(iii) u(x,0) = X for 0 < X < 1. The ‘steady-state’ or time-asymptotic limit of the solution is 
w(x,oo) = [exp{ax/n) - 1] / [exp(a/^) - Ij. 

The case a = 1, = 0.01 is selected, so that the ‘Reynolds’ number Re = a/ n = 100. 

This leads to the formation of a fairly sharp boundary layer, because the thickness of the 
layer scales as the inverse of the Reynolds number. Uniform mesh intervals Ax = 0.0025 
and At = 0.002 are used, so that the Courant number is 0.8. Fig. 9 shows the computed 
and exact steady-state limits, together with the error. The boundary layer is seen to be well 
resolved, with the maximum magnitude of the error being about 1% of the solution peak. 


8 Conclusions guid Discussions 

In the explicit a-^i scheme [1, 10], the diffusion term in Eq. (1.2) is not modeled, i.e., Eq. (2.3) 
is assumed. Also the diffusion term in Eq. (3.1) is modeled with no interpolation or extrap>- 
olation, with a resulting reduction of time-accuracy. Obviously, such a solver can be used 
only when the diffusion term is small compared with the convection term. Contrarily, the 
diffusion terms in both Eqs. (1.2) and (3.1) are modeled in the current implicit solvers. 

The current implicit solvers are carefully constructed so that they become identical to 
the explicit dual a scheme for the pure convection equation, when the viscosity coefficient 
vanishes. Because the dual a scheme is nondissipative, this construction ensures that the 
physical dissipation is never overwhelmed by numerical dissipation in the present implicit 
solvers. 

The implicit schemes have been shown to be stable provided the Courant number does 
not exceed unity in magnitude. There is no dependence of the stability of the schemes on 
the viscosity parameter, other than the condition that /i > 0. Stability analysis reveals the 
remarkable facts that (i) if // = 0, the amplification factors of the dual-mesh implicit schemes 
reduce to those of the classical Leapfrog scheme, and (ii) if a = 0, one of the amplification 
factors of the implicit a-fi schemes reduces to that of the classical Crank-Nicolson scheme. 

The truncation error analysis of the discretized equations of the implicit schemes shows 
that, if the Courant number remains bounded when refining the space-time mesh, they are 
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(;(;tisist;(;nt with the couvcjction-diffusion cc[uatioti, and arc second-order <vccurat() in space 
and time. 

Numerical examples have borne out the conclusions of the stability and truncation error 
analysis. The current implicit schemes were seen from the numerical examples to enable 
stable accurate computations over the whole viscosity range, from the pure diffusion case to 
convection dominated problems. 
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(b) an initial-value/boundary-value problem 




The space-time mesh of the implicit schemes 
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Figure 5 a-\i (II) Scheme : Sine Wave Example 
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Figure 6 Sine Wave Example : Comparison of Schemes 
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Figure 9 — a-ji (71) Scheme : Boundary Layer, Re = 100 


